Simulation of rtf = 3 QCD by Hybrid Monte Carlo 
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Simulations of odd flavors QCD can be performed in the framework of the hybrid Monte Carlo algorithm where 
the inverse of the fermion matrix is approximated by a polynomial. In this exploratory study we perform three 
flavors QCD simulations. We make a comparison of the hybrid Monte Carlo algorithm and the R-algorithm which 
also simulates odd flavors systems but has step-size errors. We find that results from our hybrid Monte Carlo 
algorithm are in agreement with those from the R-algorithm obtained at very small step-size. 



1. Introduction 

Recent lattice QCD simulations include effects 
of dynamical fermions. Due to the algorithmic 
limitation of the standard Hybrid Monte Carlo 
(HMC) algorithm Q, those simulations are lim- 
ited to even numbers of degenerate flavors. In or- 
der to include dynamical effects correctly, simula- 
tions of QCD with three flavors (u,d,s quarks) are 
desirable. Simulations with an odd number of fla- 
vors can be performed using the R-algorithm ^ . 
This algorithm, however, is not exact: it causes 
systematic errors of order At 2 , where At is the 
step-size of the Molecular Dynamics evolution. A 
careful extrapolation to zero step-size is therefore 
needed to obtain exact results. Nevertheless, it 
is common practice to forego this extrapolation 
and to perform simulations with a single step-size 
chosen small enough that the expected systematic 
errors are smaller than the statistical ones. We 
want to point out that there is an alternative to 
the R-algorithm, which gives arbitrarily accurate 
results without any extrapolation^. 

Luscher proposed a local algorithm, the so- 
called "Multiboson algorithm" Q], in which the 
inverse of the fermion matrix is approximated 
by a suitable Chebyshev polynomial. Originally 
he proposed it for two flavors QCD. Borigi and 
de Forcrand [p| noticed that the determinant 
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of a fermion matrix can be written in a mani- 
festly positive way using a polynomial approxi- 
mation, so that one can simulate odd flavors QCD 
with the multiboson method. Indeed, using this 
method, one flavor QCD was simulated success- 
fully The same polynomial approximation 
can be applied for the HMC j7|. Actually, in 
the development stage of Ref. || , one flavor QCD 
was also simulated by HMC and it was confirmed 
that the two algorithmically different methods — 
multiboson and HMC — give the same plaque- 
tte value Here we give the formulation of 
the HMC algorithm with odd flavors and perform 
rif = 3 QCD simulations. Then we compare our 
results with those of the R-algorithm. 

2. Formulation 

2.1. n f = 2 

The application of Liischer's idea to nf = 2 
QCD HMC was first made by the authors of 
Ref.@, and later by §. The lattice QCD par- 
tition function with n f = 2 degenerate quark fla- 
vors is given by 

Z = J dUdetD 2 exp(~S gauge ), (1) 

where D is the fermion matrix and in this study 
we use Wilson fermions. In the formulation of the 
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HMC algorithm, det D 2 is treated as 

det£> 2 ~ J d(/)' i d(t)exp(-<j) t D t - 1 D- 1 (t>), (2) 

where the 75 hermiticity of the fermion matrix D, 
i.e. D — 75/5^75, is used. 

Introducing momenta P conjugate to the link 
variables U , the partition function is rewritten as 

Z = J dUdPexp(-H), (3) 

where the Hamiltonian H is defined by 

H = Y -P 2 + S„ + tfD^D- 1 ^ (4) 

This Hamiltonian is used for the Molecular Dy- 
namics (MD) simulation of the standard HMC al- 
gorithm. Eq.([|) has a computational difficulty in 
MD simulations since one must solve x = D^ 1 ^ 
type equations which in general take a large 
amount of computational time for a large fermion 
matrix and/or for a small quark mass. 

Following Liischer Q , the inverse of D can be 
approximated by a polynomial: 

n 

l/D*P n (D)=l[(D-Z k ), (5) 

where Z k are the roots of the polynomial P n (D). 
We choose Z k = I — exp(i2nk/(n + 1)). 

Replacing D^ 1 in eq.(|) by P n (D) we obtain 
an approximate Hamiltonian, 

H n = X -P 2 + S„ + 4> ] ~Pn{D) ] P n {D)4>. (6) 

An advantage of using H n is that no solver calcu- 
lation is required in the MD evolution. Instead, 
one needs n multiplications by the matrix D. 

H n does introduce some systematic errors from 
the polynomial approximation. For the n/=even 
case, however, these errors are easily corrected 
by using the exact Hamiltonian of eq.(||) at the 
Metropolis step |7). This guarantees that configu- 
rations will be distributed according to the exact 
measure oc det-D 2 , for any polynomial P n (D). 

However, the domain of convergence of P n (D) 
is bounded by a circle centered at (1,0) which 
goes through the origin. If all eigenvalues of D 



fall inside this domain, P n (D) converges expo- 
nentially. Otherwise, P n (D) does not converge, 
which may happen for some exceptional configu- 
rations. Our algorithm will tend to reject these 
configurations at the Metropolis step, leading to 
extremely long autocorrelation times. This do- 
main of convergence can be changed by adopting 
another approximating polynomial. However, the 
origin must be excluded. Together with connect- 
edness and conjugate symmetry of the spectrum, 
this implies that the real negative axis is always 
excluded from the domain of convergence for any 
polynomial. Configurations with real negative 
Dirac eigenvalues will be rejected by our poly- 
nomial algorithm. 

2.2. ri/ = I 

After invention of the multiboson algorithm, 
Borigi and de Forcrand ]5| noticed that a single 
det D can be treated in a manifestly positive way 
and an nj = 1 multiboson simulation was per- 
formed to study thermodynamics of nf = 1 QCD 

B- 

As before, the inverse of the fermion matrix D, 
using a polynomial of degree 2n, is approximated 
as |,| 

In 

l/D*l[(D-Z k ), (7) 

k=l 

where Z k = 1 — exp(i 27rfc/(2n + I)). Noticing 
that the Z k s come in complex conjugate pairs, 
eq. (^j) is rewritten as 

n 

l/D^l[(D-Z k )(D-Z k ). (8) 

fc=i 

Using the 75 hermiticity of the fermion matrix, 
we find that det(D - Z k ) = det(L> - Z k )^ . Thus 
the determinant of D is written as 

det(£>) ~ det(T^(D)T n (D))-\ (9) 

where T n {D) = rifc=i(-^ — Z k ), and then we ob- 
tain 

dct(D) ~ J d(t>U(t>exp(-tfTl{D)T n (D)(l)). (10) 

The term (f>^T^(D)T n (D)(l) is manifestly positive. 
Then we may define the Hamiltonian of n/ = I 
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QCD as 
H=-P 2 



S„ 



+ tfTl{D)T n {D)4>. (11) 



With this Hamiltonian there is no difficulty to 
perform HMC algorithm. To improve efficiency 
and accuracy, one may use a polynomial of lower 
degree n during the Molecular Dynamics trajec- 
tory, and a much higher degree m 3> n for the 
Metropolis step [p|JlO|l. The domain of conver- 
gence of the approximation cq.(|To|) is the same 
as for rif = 2. Exceptional configurations for 
which eigenvalues fall outside this domain will 
likewise be rejected at the Metropolis step. A 
further difficulty is that the sampled measure 
is cx det(Tj l (D)T m ( J D))- 1 , which for exceptional 
configurations differs from the desired det D, in- 
creasingly so with m. 
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Figure 1. (top): X n versus degree n; 

(bottom) :\X n — X exact \ versus degree n. 



2.3. n f = 2 + 1 

The partition function of rtf = 2 + 1 QCD is 
given by 



Z 



dU det D 2 det D exp(—S gauge ), 



(12) 



where the notations D and D are introduced to 
distinguish the two different quark masses. Using 
cq.(g) for det D 2 and eq.@ for detD, 

det D 2 det D 
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exp(-0tL>t-i7)-i0 _ Ht ) } 
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where H T = ^T 7 \(D)T n (D)<j). We define n f 
2 + 1 Hamiltonian by 



H = \p 2 



S„ 



h 7 
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Two remarks are in order: (i) as for nt = 1, 
one could use during the MD trajectory a poly- 
nomial of lower degree than for the Metropo- 
lis step; (ii) the two bosonic fields <f) and </> 
could be replaced by a single one, with action 
^T^(D)b^- 1 b- 1 T n (D)(f). For simplicity, in this 
exploratory study we use two distinct bosonic 
fields and a single approximating polynomial. 

3. Convergence 



3.1. 



To see the rate of convergence of P n (D), we 
calculate the quantity X n = <fy P^(D)P n (D)(f). 
In the limit n — > oo, X n goes to X exact = 
^t£)t- L D- 1 <f>. First, we choose X exact = rfr) 
where 77 is a random gaussian vector. Then the 
vector 4> is set to cf> = Dr\. The accuracy of 
X n is measured by the difference between X n 
and X exact . We use a random gauge configu- 
ration. Figure l:(top) shows X n versus the de- 
gree n for different quark masses. Here the same 
77 is used for each calculation of X n . X n con- 
verges to one value as n increases, but at high 
degree n, X n diverges, which can be understood 
due to the rounding errors of our computer, where 
calculations are performed with 64-bit accuracy. 
Figure 1: (bottom) shows the accuracy of X n by 
\X n — X exact \. Exponential convergence is seen 
for each quark mass, but the rate of convergence 
is slow for small quark masses. 
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Figure 2. (top): X n versus degree n; 

(bottom):|X„ — X max \ versus degree n. 
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Figure 3. (top): Plaquette of n/ = 3 flavor 
QCD on an 8 2 x 10 x 4 lattice at = 5.0 and 
at k — 0.130 as a function of degree n; (bottom): 
Real part of Polyakov loop. 



3.2. n f = 1 

We do the same analysis as for n/ = 2, but 
for nf = 1 , the value of X exact is not known. So 
we calculate the quantity X n = (j^T^{D)T n {D)(j), 
where the vector is a gaussian random vector, 
and we use a random gauge configuration. We 
assume that X n goes to a certain value in the 
limit of n — > cxo. Figure 2: (top) shows X n as 
a function of degree n. X n seems to converges 
to a certain value when the degree n increases. 
At high degree n, X n diverges as in the case of 
n f = 2. 

To see the rate of convergence, we calculate 
\X n — X max where X max is defined by X max — 
X m , m 3> n. Due to the rounding errors, we 
can not take m very large. We take a maximum 



number m where the rounding errors still do not 
appear. Figure 2: (bottom) shows \X n — X max \ as 
a function of degree n. The dips seen in the fig- 
ure are just due to the fact that at those points 
X n = X max = X m . The convergence seems to be 
exponential, but the rate of convergence is slow 
for small quark masses as in the n/ = 2 case. 

4. Simulations 

Simulations of three flavors QCD are performed 
on an 8 2 x 10 x 4 lattice at (3 = 5.0 with k = 
0.130 and 0.160. We measure the plaquette and 
Polyakov loop varying the degree n and compare 
them with those from the R-algorithm obtained 
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with a step-size At = 0.01 Figures 3 and 

4: (top) show the plaquette as a function of n at 
k = 0.130 and 0.160, respectively. Except for 
very small n, the results from the HMC algorithm 
agree with those from the R-algorithm within sta- 
tistical errors. Results of the Polyakov loop are 
shown in Figures 3 and 4: (bottom). Except for 
a small discrepancy seen in Figure 3, the results 
from the HMC algorithm are in agreement with 
those from the R-algorithm. Note that conver- 
gence is not monotonic in n. 

5. Conclusions 

We formulated an odd-flavor HMC algorithm 
using a polynomial approximation. Simulations 
of three flavors QCD were performed. We found 
that the plaquette values are consistent with 
those from the R-algorithm at very small step- 
size. In principle the HMC algorithm is able to 
simulate any flavors of QCD, with arbitrary ac- 
curacy and without extrapolation [as long as all 
Dirac eigenvalues are not real negative] . However 
the rounding errors should be under control when 
we use a large lattice or /and small quark masses 
where one may need a polynomial of high degree 
n to achieve sufficient approximation. 

This work was partially supported by the Min- 
sitry of Education, Science, Sports and Culture, 
Grant-in- Aid, No.11740159. 
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